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In the cerebral cortex, most synapses are found in the neuropil, but 
relatively little is known about their 3-dimensional organization. Using 
an automated dual-beam electron microscope that combines focused 
ion beam milling and scanning electron microscopy, we have been 
able to obtain 10 three-dimensional samples with an average volume 
of 180 jum 3 from the neuropil of layer III of the young rat somatosen- 
sory cortex (hindlimb representation). We have used specific soft- 
ware tools to fully reconstruct 1695 synaptic junctions present in 
these samples and to accurately quantify the number of synapses per 
unit volume. These tools also allowed us to determine synapse pos- 
ition and to analyze their spatial distribution using spatial statistical 
methods. Our results indicate that the distribution of synaptic junc- 
tions in the neuropil is nearly random, only constrained by the fact 
that synapses cannot overlap in space. A theoretical model based on 
random sequential absorption, which closely reproduces the actual 
distribution of synapses, is also presented. 

Keywords: dual-beam electron microscopy, focused ion beam milling/ 
scanning electron microscopy, neocortex, spatial distribution of synapses, 
spatial statistics, synapses 



Introduction 

In the cerebral cortex, most synapses (90-98%) are established 
in the neuropil (Alonso-Nanclares et al. 2008), which is made 
of dendrites, axons, and glial processes. One major issue in 
cortical circuitry is to ascertain how synapses are distributed 
and whether synaptic connections are specific or not (DeFe- 
lipe et al. 2002). Therefore, to understand the anatomical 
design principles of the cortical circuits, it is essential to 
analyze the ultrastructure of all components of the neuropil 
and, in particular, the number and spatial distribution of sy- 
napses. There are 2 major morphological types of synapses: 
Asymmetric (or Gray's type I) and symmetric (or Gray's type 
II; Gray 1959; Colonnier 1968). The major sources of asym- 
metric synapses are spiny neurons (pyramidal and spiny non- 
pyramidal cells) and extrinsic cortical afferents, whereas the 
vast majority of symmetric synapses are of intrinsic origin and 
are formed by the population of aspiny or sparsely spiny inter- 
neurons. Since spiny neurons and the major cortical afferent 
systems are excitatory in function, and virtually all aspiny or 
sparsely spiny interneurons are y— aminobutyric (GABA acid)- 
ergic, it is assumed that axon terminals forming asymmetric 
synapses are excitatory and those forming symmetric synapses 
are inhibitory (Ascoli et al. 2008). Furthermore, synaptic size 



plays an important role in the functional properties of sy- 
napses (Schikorski and Stevens 1997; Takumi et al. 1999; 
Liischer et al. 2000; Tarusawa et al. 2009)- Thus, numerous re- 
searchers have been trying to find simple and accurate 
methods for estimating the distribution, size, and number of 
synapses. To this end, 2 sampling procedures are currently 
available: One is based on serial reconstructions and the other 
on single sections. Clearly, serial reconstruction should be the 
method of choice for the challenging task of unraveling the 
extraordinary complexity of the nervous system. Indeed, serial 
sectioning transmission electron microscopy is a well estab- 
lished and mature technology for obtaining 3-dimensional 
data from ultrathin sections of brain tissue (Stevens et al. 
1980; Harris et al. 2006; Hoffpauir et al. 2007; Mishchenko 
et al. 2010; Bock et al. 2011). It is based on imaging ribbons 
of consecutive sections with a conventional transmission 
electron microscope. However, the major limitation is that ob- 
taining long series of ultrathin sections is extremely time- 
consuming and difficult, often making it impossible to recon- 
struct large volumes of tissue. Hence, the recent development 
of automated electron microscopy techniques represents 
another crucial step in the study of neuronal circuits (e.g. 
Briggman and Denk 2006; Knott et al. 2008; Merchan-Perez, 
Rodriguez, Alonso-Nanclares, et al. 2009)- 

In the present study, we use a new dual-beam electron 
microscope that combines a focused ion beam (FIB) column 
and a scanning electron microscope (SEM). The FIB column 
directs a gallium (Ga + ) ion beam toward the sample. The col- 
lision of Ga + ions with substrate atoms results in the removal 
of these atoms from the substrate. Since the FIB can be 
focused and controlled on a nanometer scale, this effect can 
be used to mill the specimen, removing a thin layer of 
material. A backscattered electron image is then obtained from 
the freshly milled surface using the SEM. The milling/imaging 
processes are automatically repeated in order to obtain a 
long series of images that represent a 3-dimensional sample of 
tissue. Using this methodology, we have previously shown 
that virtually all cortical synapses can be accurately identified 
as asymmetric (excitatory) or symmetric (inhibitory) synaptic 
junctions when they are visualized, reconstructed, and quanti- 
fied from large 3-dimensional tissue samples obtained in an 
automated manner (Merchan-Perez, Rodriguez, Alonso- 
Nanclares, et al. 2009)- Since the segmentation, quantification, 
and analysis of synaptic junctions are a labor-intensive pro- 
cedure, we used Espina software, a specifically developed tool 
that greatly facilitates and accelerates these processes (Morales 
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et al. 2011). For the segmentation of synaptic junctions, 
Espina makes use of the fact that presynaptic and especially 
postsynaptic densities are electron-dense structures. It uses a 
gray-level threshold to extract all the voxels that fit the gray 
levels of the synaptic junction. Asymmetric synaptic junctions 
have a thicker postsynaptic density, so they can be distin- 
guished from symmetric ones, whose postsynaptic densities 
are thinner and similar to the presynaptic densities. As a 
result, each synaptic junction (already tagged as asymmetric 
or symmetric) is segmented as a flattened irregular volume 
and its geometrical features, including its spatial position, are 
determined by the same software. 

Points in space — the synaptic junctions in the present case — 
can be distributed according to 3 basic patterns: Random, 
clustered, and regular. A random pattern follows the basic re- 
ference model of a point process, the so-called complete 
spatial randomness (CSR) or homogeneous spatial Poisson 
point process (Seidel 1876). It is characterized by 2 funda- 
mental properties: (1) The number of points falling in a given 
region follow a Poisson distribution with a constant mean 
number of points per unit volume; and (2) each point is 
equally likely to occur at any position and is not affected by 
the locations of the other points (Illian et al. 2008; Gaetan 
and Guyon 2009). In a clustered distribution, the points tend 
to concentrate in groups, leaving other regions of space 
empty or almost devoid of points. Finally, in a regular or dis- 
persed pattern, every point is located as far as possible from 
its neighbors, resulting in a lattice-like distribution of points. 
Although there are no clear-cut limits between these 3 basic 
patterns, CSR represents a boundary condition between clus- 
tered and dispersed spatial processes, so our strategy was first 
to test for CSR, trying to identify deviations from it and to 
build a model that accounted for these deviations (Illian et al. 
2008; Gaetan and Guyon 2009). 

Materials and Methods 

Tissue Preparation and 3-Dimensional Electron Microscopy 

Male Wistar rats sacrificed in postnatal day 14 were used for this 
study. Animals were administered a lethal intraperitoneal injection of 
sodium pentobarbital (40 mg/kg) and were intracardially perfused 
with 2% paraformaldehyde and 2.5% glutaraldehyde in 0.1 M phos- 
phate buffer. The brain was then extracted from the skull and pro- 
cessed for electron microscopy according to a previously described 
protocol (Merchan-Perez, Rodriguez, Alonso-Nanclares, et al. 2009). 
All animals were handled in accordance with the guidelines for 
animal research set out in the European Community Directive 2010/ 
63/EU, and all procedures were approved by the local ethics commit- 
tee of the Spanish National Research Council (CSIC). 

Three-dimensional brain tissue samples were obtained using a com- 
bined FIB/SEM (Crossbeam® Neon40 EsB, Carl Zeiss NTS GmbH, 
Oberkochen, Germany). This instrument combines a high-resolution 
field emission SEM column with a focused gallium ion beam, which 
can mill the sample surface, removing thin layers of material on a nan- 
ometer scale. After removing each slice (20 nm of thickness), the 
milling process was paused, and the freshly exposed surface was 
imaged with a 1.8-kV acceleration potential using the in-column 
energy selective backscattered (EsB) electron detector. A 30-mm aper- 
ture was used and the retarding potential of the EsB grid was 1500 V. 
The milling and imaging processes were sequentially repeated, and 
long series of images were acquired through a fully automated pro- 
cedure, thus obtaining a stack of images that represented a 
3-dimensional sample of the tissue (Merchan-Perez, Rodriguez, 
Alonso-Nanclares, et al. 2009). Image resolution in the xy plane was 
3 7 nm/pixel. Resolution in the z axis (section thickness) was 20 nm. 



Ten different samples of the neuropil of layer III of the somatosen- 
sory cortex were obtained from 3 different animals (Table 1). To select 
the exact location of the samples, we first obtained plastic semithin 
sections (2 pm thick) from the block surface and stained them with 
toluidine blue to identify cortical layers. These sections were then 
photographed with a light microscope. The last of these light micro- 
scope images (corresponding to the section immediately adjacent to 
the block face) was then collated with SEM photographs of the surface 
of the block. In this way, it was possible to accurately identify the 
regions of the neuropil to be studied. All samples selected for FIB 
milling/SEM imaging were located at mid-depth of layer III. After ap- 
plying a 3-dimensional unbiased counting frame and correcting for 
tissue shrinkage (Howard and Reed 2005; Merchan-Perez, Rodriguez, 
Alonso-Nanclares, et al. 2009; Morales et al. 2011), the volume of 
tissue analyzed from each sample ranged from 149.13 to 247.58 pm 3 
(mean 180.33 pm 3 ; see Supplementary Table 1). Synaptic junctions 
within these volumes were visualized, automatically segmented, and 
reconstructed in 3 dimensions with Espina software (Morales et al. 
2011; Fig. 1). This software also provided data such as the volume of 
the unbiased counting frame, the number of synaptic junctions inside 
the frame, the spatial position of the centroids or centers of gravity of 
the synaptic junctions, and an estimation of their sizes by the Feret's 
diameter (the diameter of the smallest sphere circumscribing the sy- 
naptic junction; Fig. 1 and Supplementary Videos 1 and 2). As pre- 
viously described in detail (Merchan-Perez, Rodriguez, Alonso- 
Nanclares, et al. 2009), synaptic junctions with either a prominent or 
thin postsynaptic density were classified as asymmetric or symmetric 
synaptic junctions, respectively (Supplementary Video 1). 



Spatial Point Pattern Analysis 

To analyze the spatial distribution of synaptic junctions, we compared 
the actual positions of the centroids of synaptic junctions with 2 
theoretical point processes: CSR and random sequential adsorption 
(RSA). CSR defines a situation where a point is equally likely to occur 
at any location within the study volume, regardless of the locations of 
other points (Diggle 2003). The mathematical model underlying CSR 
is a homogeneous spatial Poisson point process. Thus, the number of 
points occurring within a finite volume, V, is a random variable fol- 
lowing a Poisson distribution with mean X \ V\ , where A > 0 and | V\ is 
the volume of V. The intensity A represents the expected number of 
points per unit volume. CSR serves as a boundary condition between 
spatial processes that are more clustered than random and spatial pro- 
cesses that are more regular than random. 

For modeling the spatial distribution of synapses, we should take 
into account that synaptic junctions cannot overlap, and thus the 
minimum intersynaptic distances (measured between their centroids) 
must be limited by the sizes of synaptic junctions themselves. A 



Table 1 

Data on the spatial distribution of synaptic junctions in 10 samples of the neuropil of layer III 



Sample no. and animal 


Number of 


Mean distance to 


Mean Feret's diameter 


identification 


synapses per 


nearest neighbor 


of synaptic junctions 




cubic micron 


(nm) ± SD 


(nm) ± SD 


1 (W31) 


0.9857 


519.55 ± 136.35 


377.19 ± 159.63 


2 (W31) 


0.6936 


594.07 ± 192.28 


462.18 ± 177.52 


3 (W33) 


0.9279 


537.43 ± 159.20 


437.62 ± 168.04 


4 (W33) 


1.0088 


537.39 ± 157.70 


414.22 ± 169.04 


5 (W33) 


0.9474 


597.30 ± 174.02 


466.03 ± 215.91 


6 (W33) 


0.9399 


533.21 ± 163.29 


423.38 ± 169.83 


7 (W33) 


0.9881 


487.17 ± 172.30 


397.29 ± 168.22 


8 (W35) 


0.7997 


568.21 ± 178.51 


427.90 ± 168.15 


9 (W35) 


1.1267 


501.38 ± 156.97 


378.35 ± 166.60 


10 |W35) 


1.0178 


523.74 ± 150.36 


405.43 ± 175.62 


All samples 


0.9399 


535.78 ± 166.81 


417.06 ± 175.97 



Note: Densities of synapses were calculated by dividing the actual number of synaptic junctions 
within an unbiased 3-dimensional counting frame by the volume of the counting frame. Distances 
to the nearest neighbor were calculated between the centroids of the synaptic junctions. 
The Feret's diameter is the diameter of the smallest sphere circumscribing the synaptic junction; 
it is used here as an estimate of the size of synaptic junctions. 
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Figure 1. Segmentation of synaptic junctions, estimation of their sizes, and extraction of their geometric centers or centroids. (4) An example of the tissue volume sampled 
and the reconstructed synaptic junctions. Asymmetric and symmetric synaptic junctions are shown in green and red, respectively. (fl) The smallest spheres circumscribing the 
reconstructed synaptic junctions served to calculate the Feret's diameter, which was used as an estimation of their size. An unbiased counting frame was also drawn to facilitate 
the quantification of the number of synapses per unit volume. (C) The geometric centers or centroids of the synaptic junctions were determined to indicate the spatial position of 
the synapses. They are represented within the spheres that were used to estimate the sizes of synaptic junctions. (0) The centroids indicating the spatial position of synaptic 
junctions are represented alone within the counting frame. The dimensions of the counting frame represented in B, C, and D are 7.16 x 4.58 x 3.98 /jm. See also 
Supplementary Videos 1 and 2. 
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Figure 2. Representation of the distribution of the sizes of synaptic junctions 
estimated by their Feret's diameters in the 10 samples of the neuropil of layer III 
(blue histogram) and the corresponding log-normal distribution (red trace). 

spatial point process to model this situation can be seen as a type of 
hard-core point process where no points can occur at a distance 
smaller than a specific minimum distance (Mian et al. 2008). In this 



respect, an RSA process (Evans 1993) is a kind of hard-core process 
where the pattern is constructed by placing spheres in 3-dimensional 
space iteratively and randomly, with their radii following a probability 
density function. If any newly generated sphere intersects with an ex- 
isting sphere, the new sphere is rejected and another sphere with a 
different center and radius is generated. This process is stopped when 
the required number of spheres is reached (Supplementary Video 3). 

In our case, the minimum interpoint distance would depend at 
least partially on the sizes of synaptic junctions. Thus, we obtained the 
empirical distribution of the Feret's diameters of synaptic junctions as 
an estimate of their sizes. We then performed a goodness-of-fit test to 
find the theoretical probability density distribution that best fitted to 
the empirical distribution. We found that the best fit using the Kolmo- 
gorov-Smirnov test was with a log-normal distribution with 
parameters fi = 5.828 and (7 = 0.446 (Fig. 2). The estimates of the 
expectation and variance were v = 375.198 nm and 1^ = 30,981.351 
nm 2 , respectively (Supplementary Methods). 

For each of the 10 different samples and their corresponding CSR 
and RSA processes, we calculated 3 functions commonly used for 
spatial point pattern analysis: the G, F, and K functions (Ripley and 
Kelly 1977; Diggle 1979). 

The G function, also called the nearest-neighbor distance cumulat- 
ive distribution function or the event-to-event distribution, is, for a 
distance d, the probability that a typical point separates from its 
nearest-neighbor a distance of at most d. The empirical distribution 
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function G* of the observed nearest-neighbor distances is used to esti- 
mate the G function: 

where D t is the distance of synapse i to its nearest-neighbor synapse, 
N the number of synapses, and /(■) the indicator function which is 
one if the argument is true and is zero otherwise. 

The F function, also called the empty space function or the 
point-to-event distribution, is, for a distance d, the probability that the 
distance of each point (in a regularly spaced grid of L points super- 
imposed over the sample) to its nearest synapse centroid is at most d. 
The empirical distribution function F* of the observed empty space 
distances on a grid of L locations is used to estimate the F function: 

F*(d)=i^/(£>;<rf), 

i 

where D\ is the distance of grid point i to its nearest centroid. No 
edge correction was used for G and F functions. 

The K function, also called the reduced second moment function 
or Ripley's function, is, for a distance d, the expected number of 
points within a distance d of a typical point of the process divided by 
the intensity X. An estimation of the K function is given by the mean 
number of points within a sphere of increasing radius d centered on 
each sample point, divided by an estimation of the expected number 
of points per unit volume. A common formula is the Milles-Lantue- 
joul-Stoyan-Hanisch translation-corrected estimate (Baddeley et al. 
1993), given by: 

where | V\ is the volume of the counting frame, /?"« is the distance 
between centroid i and centroid j, and | Vg | denotes the volume of 
the 3-dimensional subspace obtained as the intersection between V 
and the (translated) sphere centered at synapse i and radius D" y. 

Spatial analysis of the positions of the centroids of synaptic junc- 
tions in experimental and simulated distributions was performed with 
R Project, specifically with "spatstat" (Baddeley and Turner 2005) and 
"fitdistplus" (Venables and Ripley 2002; Vose 2008) packages and 
with SA3D software (Eglen et al. 2008). Additional statistical analyses 
were performed with SPSS (IBM Corp., New York, USA). 



Results 

We obtained 10 independent samples by FIB/SEM technology 
from regions of the neuropil of layer III of the somatosensory 
cortex (hindlimb representation) from three 14-day-old rats. 
From each of the 10 samples of cortical tissue used in this 
study, we obtained the number, position, and shape of synap- 
tic junctions inside a 3-dimensional unbiased counting frame 
(Howard and Reed 2005; Fig. 1 and Supplementary Video 2). 
No cell somata or blood vessels were present within the count- 
ing frames. In this way, we directly calculated the number of 
synapses per unit volume of neuropil. These densities ranged 
from 0.69 to 1.13 synapses/um 3 , with a mean of 0.94 sy- 
napses/um 3 (Table 1). The spatial positions of the geometric 
barycenters or centroids of each individual synaptic junction 
(Fig. 1 and Supplementary Video 2) were recorded for further 
analysis. As an estimation of the size of synaptic junctions, we 
also recorded the Feret's diameter of each synaptic junction, 
which is the diameter of the smallest sphere circumscribing 
the synaptic junction (Fig. 1, Table 1, and Supplementary 
Videos 1 and 2). The total tissue volume studied for the 10 



samples was 1803-32 um 3 , where we found 1695 synapses 
(1555 asymmetric and 140 symmetric). Additional data on 
these tissue volumes and the actual counts of asymmetric and 
symmetric synaptic junctions per sample are shown in Sup- 
plementary Table 1 . 

In order to explore the spatial distribution pattern of synap- 
tic junctions, we first performed a nearest-neighbor analysis 
on the centroids of the synaptic junctions found within each 
sampled volume. We calculated the mean distance from each 
centroid to its nearest neighbor within the counting frame. 
Centroids that were closer to the boundaries of the counting 
frame than to any other centroid were excluded from the cal- 
culations, since their nearest neighbor could lay outside the 
counting frame at an unknown distance (Baddeley et al. 1993; 
Illian et al. 2008). The mean distance to the nearest neighbor 
ranged from 487.17 to 597.30 nm (Table 1). 

To further analyze the spatial distribution of nearest neigh- 
bors, we obtained the G function (Diggle 1979). This function 
is estimated using the distances from each centroid to its 
nearest neighbor, and plotting the fraction of points in the 
sample having their nearest neighbor at a given distance or 
less (Fig. 3). We then compared the experimentally observed 
G functions with the theoretical G functions that would be ob- 
tained if the same number of centroids were randomly distrib- 
uted according to the CSR or homogeneous Poisson process. 
In the latter conditions, the function takes a sigmoid shape, as 
in the examples shown in Figure 3A,C (red traces) and Sup- 
plementary Figures 1-10. In a clustered pattern, the curve 
would be skewed to the left (shorter interpoint distances than 
expected), while it would be shifted to the right in a regular 
pattern (longer interpoint distances than expected). 

The experimentally observed G functions showed a similar 
shape in the 10 different samples. Two examples of them are 
shown in Figure 3 (blue traces) and the rest are shown in Sup- 
plementary Figures 1-10. In all cases, when the experimen- 
tally observed G functions were compared with the G 
functions of the corresponding Poisson processes, the more 
noticeable differences were found at short interpoint dis- 
tances, since all experimental samples showed an empty space 
around centroids where no other centroids were found 
(arrows in Fig. 3^4, C, see also Supplementary Figs 1-10). Due 
to this empty space, the fraction of centroids having their 
nearest neighbor at short distances was lower than would be 
expected in a homogeneous Poisson process, while at longer 
distances the observed curves tended to overlap or slightly 
exceed the Poisson curves. In general, 20-40% of the centroids 
had their nearest neighbor at 500 nm or less, while practically 
all of them had their nearest neighbor at 1000 nm or less 
(Fig. 3A, C and Supplementary Figs 1-10). 

The existence of this dead space can be explained by the 
fact that the centroids of synaptic junctions represent the geo- 
metric center of a physical volume (the volume of the recon- 
structed synaptic junction) instead of a dimensionless point, 
and thus their position cannot be independent from the pos- 
ition of their immediate neighbors. In other words, centroids 
cannot be too close to each other since the volumes they rep- 
resent, the synaptic junctions, cannot overlap. Therefore, in 
order to simulate the distribution of synapses in the neuropil, 
the empty space around the centroids of synaptic junctions 
had to be accounted for by a random minimum-spacing 
distance that would depend, at least partially, on the sizes of 
synaptic junctions themselves. 
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Figure 3. G functions estimated from the positions of the centroids of synaptic junctions of 2 examples of the neuropil of layer III, corresponding to samples no. 5 ffi and B) and 
no. 10 (C and D) (see also Table 1). \nA and C, the functions obtained from the experimentally observed data (blue traces) have been represented together with the functions 
that would be obtained if the centroids were distributed according to a homogenous Poisson process (red traces) and the mean of 100 simulated realizations of an RSA process 
(green traces, see text for further details). In both samples, the experimentally observed functions show an empty space around centroids where no other centroids are found 
(4 and B, arrows). A similar empty space is also present in RSA simulations. In B and D, the experimentally observed G functions are depicted in blue and the 100 individual 
simulated realizations of RSA in green, showing that the experimental G functions either overlapped with the simulated RSA functions (6) or were slightly displaced to the left 
(D). See also Supplementary Figures 1-10. 



To account for this minimum-spacing distance between cen- 
troids, we simulated random distributions of points in space 
using a sequential algorithm that first generates a random 
point that is the center of a sphere with a certain diameter. The 
diameters used in the simulations were randomly drawn from 
the probability distribution of the experimentally observed 
Feret's diameters (Table 1). This distribution was found to be 
log-normal (Fig. 2). If the next simulated point with its corre- 
sponding Feret's diameter does not overlap with previously 
generated points, then it is accepted. Otherwise, it is rejected 
and another random point is generated. The process termi- 
nates when the desired number of points has been reached 
(Supplementary Video 3). This kind of process is known as 
RSA (Evans 1993). 

In order to compare the experimentally observed distribution 
of centroids with the RSA process, 100 simulated realizations of 
RSA were generated for each experimental sample, with the 
same number of points and volumes as the original samples. 
Fitness between the observed G functions and the ones 
obtained from simulated RSA processes was better than that 
previously found with homogeneous Poisson processes, 
although the observed curves tended to be located slightly to 
the left of the RSA curves (Fig. 3 and Supplementary Figs 1-10). 

Further analysis of the spatial distribution of the centroids 
of synaptic junctions was performed by calculating the F and 
K functions (Baddeley et al. 1993; O'Sullivan and Unwin 



2002; Gaetan and Guyon 2009). To estimate the F function, a 
regular grid is traced within the 3-dimensional bounding box 
that contains the centroids, the distances between each grid 
crossing point and its nearest neighboring centroid are 
measured, and the cumulative probability of having the 
nearest centroid at a given distance or less is plotted. In a 
homogeneous Poisson process, the F function takes the same 
values as the G function (Baddeley et al. 1993; Fig. 4 and Sup- 
plementary Figs 1-10). Values of the F function greater than 
the Poisson value suggest that there is regularity in the point 
pattern, while lower values suggest clustering. The K function 
is estimated as the mean number of points within a sphere of 
increasing radius centered on each sample point (Ripley and 
Kelly 1977). It is a rapidly growing curve in a homogeneous 
Poisson process (Fig. 5 and Supplementary Figs 1-10), with 
lower values in a regular pattern and higher values in a clus- 
tered pattern. 

In the case of F functions, the Poisson and RSA curves 
were very similar, and the experimentally observed curves 
either overlapped with RSA curves or tended to be located 
slightly to the right of them (Fig. 4 and Supplementary Figs 
1-10). K functions were also very similar for both the Poisson 
and RSA models. The experimentally observed curves of the 
K function were very similar to the Poisson and RSA curves or 
were located slightly to the left of them (Fig. 5 and Sup- 
plementary Figs 1-10). 
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Figure 4. F functions corresponding to the same samples (5 and 10) as in Figure 3. Colors of traces represent the same as in Figure 3. The F functions were similar for 
the experimentally observed data and the Poisson and RSA processes. The observed functions either overlapped with ifl) or were slightly displaced to the right (C) of the 
corresponding Poisson and RSA curves. Note that in both examples the observed functions were compatible with individual realizations of the RSA model (fl and D). 
The F functions of the rest of the samples are shown in Supplementary Figures 1-10. 



In order to test the goodness-of-fit of the experimentally 
observed functions in the 10 different samples with those for 
the corresponding Poisson and RSA processes, we performed 
a battery of 2-sample Kolmogorov-Smirnov tests (Table 2). 
The tests showed that fitness was better for F and K functions 
than for G functions. Comparing the actual samples with RSA 
simulations yielded a clearly better fit than comparing with 
Poisson distributions. In fact, in the case of RSA simulations, 
26 of 30 tests yielded P-values higher than 0.950, while for 
Poisson distributions 17 of 30 tests yielded P-values over 
0.950 (Table 2). Only for sample 2, the null hypothesis of RSA 
was rejected for the G function (P = 0.035). For the homo- 
geneous Poisson process, the null hypothesis was only re- 
jected for the G function of sample 10 (P= 0.007). In these 
samples, however, the G functions showed no clear deviations 
toward either a clustered or regular distribution pattern (Sup- 
plementary Figs 2 and 10). 

Results were similar when all synapses (asymmetric and 
symmetric) were studied as a single group and when only 
asymmetric synapses were analyzed (Supplementary Table 2). 
This result was expected since asymmetric synapses rep- 
resented 91-74% of the total number of synapses found in 
layer III. However, due to the small number of symmetric sy- 
napses present in our samples (the average number of sym- 
metric synaptic junctions per sample was 14, see 
Supplementary Table 1), it was not possible to calculate the 
G, F, and K functions for them. 

There is a simple way, however, to test whether asymmetric 
and symmetric synapses are intermingled at random. If they 



were intermingled at random (null hypothesis), the nearest 
neighbor of any given synapse should be asymmetric or sym- 
metric with a probability that would only depend on the 
respective proportions of both types of synapses in the 
sample. On the other hand, if both types of synapses were 
spatially segregated, their respective nearest neighbors would 
tend to be of the same type. Note that in this case we did not 
measure the distances between centroids but counted how 
many synapses had a nearest neighbor of the same or differ- 
ent type. These counts showed that 704 asymmetric synapses 
had an asymmetric nearest neighbor and 38 had a symmetric 
nearest neighbor. In the case of symmetric synapses, 43 had 
an asymmetric nearest neighbor, while 20 had a symmetric 
one. These data were used to create a 2 x 2 contingency table 
showing both types of synapses against the type of their 
nearest neighbor (Table 3). The 2-tailed Fisher's exact test 
applied to the contingency table rejected the hypothesis that 
both types of synapses are intermingled at random 
(P< 0.0001). In fact, the number of symmetric synapses that 
have a symmetric nearest neighbor was more than 4 times 
higher than expected under the null hypothesis (Table 3), 
possibly indicating that symmetric synapses tend to group 
together. 



Discussion 

Using FIB/SEM technology, we have been able to accurately 
quantify the density and spatial distribution of synapses in 
stacks of images from the neuropil of layer III of the young rat 
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Figure 5. /(functions corresponding to the same samples (5 and 10) as in Figures 3 and 4. Trace colors represent the same as in Figure 3. The experimentally observed curves 
of the K function were generally located slightly to the left of the simulated RSA curves, although the observed functions were compatible with individual realizations of the RSA 
process. The remaining samples are shown in Supplementary Figures 1-10. 



Table 2 

P-values of the 2-sample Kolmogorov-Smirnov tests performed with the G, F and K functions 
derived from the actual experimental samples, from homogeneous Poisson processes, and from 
the mean of 100 simulated realizations of RSA per sample 



Sample no. Two-sample Kolmogorov-Smirnov tests (P-values) 

Experimental samples versus Experimental samples versus 

homogeneous Poisson processes simulated realizations of RSA 





G 


F 


K 


G 


F 


K 


1 


0.952 


0.952 


0.468 


0.799 


0.952 


0.994 


2 


0.799 


0.999 


0.699 


0.035 


0.999 


0.906 


3 


0.998 


1.000 


0.994 


0.952 


0.999 


0.994 


1 


0.998 


0.999 


0.906 


0.998 


1.000 


0.994 


5 


0.998 


1.000 


0.699 


0.999 


1.000 


0.994 


6 


0.952 


0.999 


0.906 


0.952 


0.999 


0.994 


7 


0.236 


0.952 


0.994 


0.134 


0.998 


0.994 


8 


0.998 


0.999 


0.906 


0.998 


1.000 


0.999 


9 


0.071 


0.134 


0.994 


0.998 


0.998 


1.000 


10 


0.007 


0.236 


0.906 


0.998 


0.952 


1.000 



somatosensory cortex. In addition, virtually all synapses 
present in the 3-dimensional samples could be unambigu- 
ously identified since they can be visualized in consecutive 
serial sections and, if necessary, they can be digitally resec- 
tioned at different planes to ascertain their identity as 
asymmetric or symmetric (Merchan-Perez, Rodriguez, /Alonso- 
Nanclares, et al. 2009). In this way, the technology used in the 
present study eliminates the relatively large pool of synapses 
that were labeled as "uncharacterized" in previous reports 
since they could not be classified as asymmetric or symmetric 



Table 3 

Contingency table showing the type of synaptic junctions against the type of their nearest 
neighbors 



Type of synapse 


Type of nearest neighbor 






Asymmetric 


Symmetric 


Total 


Asymmetric 


704 


38 


742 




688.54 


53.46 


742.00 


Symmetric 


43 


20 


63 




58.46 


4.54 


63.00 


To la I 


747 


58 


805 




747.00 


58.00 


805.00 



The observed counts of synapses (bold) were obtained from the 1 0 samples of the neuropil of 
layer III. The expected counts (in italics) are calculated from the marginal totals. Synapses whose 
nearest neighbors were not known — those that were located closer to the boundaries than to 
any other synapse in the sample — were excluded from the analysis. The Fisher's exact test 
applied to the observed counts rejected the hypothesis that asymmetric and symmetric synapses 
are intermingled at random (P < 0.0001). The most salient mismatch between the observed and 
expected values was the number of symmetric synapses that had a symmetric nearest neighbor, 
which was 4.4 times higher than expected. 



from single 2-dimensional images (see for example DeFelipe 
et al. 1999). Unfortunately, the mean density of synapses per 
unit volume that we have obtained in this study (mean of 0.94 
synaptic junctions/um 3 ) and the proportion of asymmetric 
and symmetric synapses could not be compared with previous 
studies since no detailed electron microscope data are avail- 
able in this particular layer of the rat somatosensory cortex at 
the age analyzed in the present study. 
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Regarding the spatial position of synapses, our results indi- 
cate that the distribution of synaptic junctions in the neuropil 
is nearly random, only constrained by the fact that synapses 
cannot overlap in space and so their geometric centers or cen- 
troids cannot be too close to their neighbors. This conclusion 
is further supported by the finding that the spatial distribution 
of the centroids of synaptic junctions can be closely modeled 
by a random process where overlapping is prevented by as- 
signing a volume to each centroid. In mathematical terms, this 
model is known as "random sequential adsorption" (Evans 
1993; Stoyan and Schlather 2000; van Lieshout 2006). None of 
the functions F, G, and K alone suffice for the characterization 
of a point pattern, so we have used them in combination 
(O'Sullivan and Unwin 2002; Gaetan and Guyon 2009). Visual 
inspection of all these functions showed only small deviations 
between the experimental data and the RSA model, which was 
further corroborated by a battery of goodness-of-fit tests. In 
our particular case, the volumes that were randomly assigned 
to each centroid were not arbitrary, but were drawn from the 
distribution of Feret's diameters of the 3-dimensionally recon- 
structed synaptic junctions. Although the Feret's diameter is 
only a rough estimate of the sizes of synaptic junctions, the 
RSA model not only prevents the overlapping of synaptic junc- 
tions, but also very closely mimics the actual distribution of 
synapses in the neuropil. Moreover, since no other constraints 
were imposed, the positions of synapses are otherwise nearly 
independent of the position of their neighbors. This and the 
fact that the Feret's spheres occupy only 5.85% of the total 
volume of the neuropil help to explain why the distributions 
of centroids in the samples and the simulated RSA processes 
are still very similar to a CSR process. 

It must be stressed that our claim that the spatial distri- 
bution of synapses follows an RSA process does not necessary 
imply that during development synapses are actually formed 
and positioned in space following the same rules as the RSA. 
However, it is interesting to consider the idea that a synapse 
could be formed anywhere in space where an axon terminal 
and a dendritic element may touch, provided this particular 
spot is not already occupied by a preexisting synapse. Indeed, 
a similar model based on random axonal outgrowth and com- 
petition for available space has already been proposed 
(Kaiser et al. 2009). Nevertheless, the same end result could 
be achieved in different ways, including the synaptic turnover 
that is known to occur during cortical development (Rakic 
et al. 1994), provided that synapses are randomly added and/ 
or withdrawn from the population. 

Since there seems to be no limitation to the position of any 
synapse except the space already occupied by other synapses, 
at least a subpopulation of them would be located very close 
to one another. That is the case in our samples, since G func- 
tions show that approximately 20-40% of synaptic junctions 
have a nearest neighbor with an intercentroid distance of 
<500 nm. Considering the actual size of synapses (Table 1), 
this would mean that some of them would be located side by 
side, opening the possibility that the neurotransmitter released 
by one synapse could reach its neighbor by diffusion and 
influence its behavior (Fuxe et al. 2007; Rusakov et al. 2011). 
However, this seems to be true mostly for glutamatergic sy- 
napses (asymmetric synapses, the most abundant type), 
because the number of GABAergic synapses (symmetric sy- 
napses) is so low that very few of them were found near 
another asymmetric synapse. This is an interesting difference 



between the neuropil and perisomatic innervation regarding 
the organization of synaptic circuits. Indeed, there are numer- 
ous glutamatergic and perisomatic GABAergic axon terminals 
in close proximity around the soma and proximal dendrites of 
pyramidal cells, providing a putative basis for nonsynaptic 
interactions between these 2 types of synapses (Merchan- 
Perez, Rodriguez, Ribak, et al. 2009). 

An issue that needs to be addressed is whether the den- 
sities of synaptic junctions and their nearly random distri- 
bution in space have any predictive value to calculate the 
possible number of synaptic inputs that a given neuron may 
receive in layer III. This is related to the so called Peter's rule 
(Braitenberg and Schiiz 1998), which represents a generaliz- 
ation to all types of cortical synapses based on the conclusion 
of Peters and Feldman (1976) regarding the synaptic connec- 
tions of geniculo-cortical afferents and cortical neurons in the 
rat visual cortex: "Although thalamocortical afferents termi- 
nate principally in layer IV, their distribution with respect to 
postsynaptic targets may be essentially random, in the sense 
that no specific types of neurons receive the afferents. 
Instead, all elements in layer IV that are capable of forming 
asymmetric synapses may be involved" (Peters and Feldman 
1976). 

Our results seem to support this assertion since the pos- 
itions of synapses do not show any spatial preference, so the 
rule would apply to the neuropil of layer III as a whole. Our 
findings are also in line with the distribution of swellings 
along pyramidal cell axon collaterals stained by the Golgi 
method in the mouse neocortex. In that case, the intervals 
between the axon swellings (characterized as putative sy- 
napses) were found to follow a Poisson distribution, with a 
dead space caused by the spatial extent of the swellings 
(Hellwig et al. 1994). However, it must be pointed out that 
the densities of synapses do vary locally (Table 1), and that 
they are randomly distributed, so local deviations from the 
rule will arise simply by chance. In any case, it is clear that 
the formulation of Peter's rule should take into account that 
different neuron types have different bouton distributions, 
project preferentially to different types of neurons and/or to 
different parts of these neurons (dendritic tree, soma, axon 
initial segment), and that both the axon terminals and dendri- 
tic trees may or may not originate in the same cortical layer 
(Anderson et al. 2002; Binzegger et al. 2004; Mishchenko 
et al. 2010; Bock et al. 2011). Thus, our data represent a het- 
erogeneous population of synapses, and our conclusion that 
they are randomly distributed in space should be maintained 
only for this population as a whole, with the only possible 
exception of inhibitory synapses (see below). 

Spatial randomness does not necessarily mean unspecific 
connections, since this does not preclude the existence of sub- 
jacent specificities based on molecular affinities, electrical 
activity, the developmental history of synapses, or other 
mechanisms not detected by our methodology. Indeed, our 
technique is capable of distinguishing between excitatory and 
inhibitory synapses, but there is no indication as to their par- 
ticular neuronal type of origin. Different groups of the sy- 
napses in any given sample could come from cells located in 
the immediate vicinity, from other layers, from cortical areas 
outside the somatosensory cortex, from the contralateral 
cortex, or from subcortical regions. If all these different types 
of synapses were to form segregated groups (e.g. DeFelipe 
and Jones 1991), it would not be detected by our present 
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methodology. However, we did find evidence that symmetric 
synapses tend to be spatially segregated from asymmetric sy- 
napses. Although we could not apply the same tests to asym- 
metric and symmetric synapses due to the scarcity of the 
latter, we found that symmetric synapses have a nearest 
neighbor of the same type over 4 times more than would be 
expected if they were intermingled at random with asym- 
metric synapses. A possible explanation is that these synapses 
may originate from one or few single axons, giving "en 
passant" synapses or ramifying locally, as typically occurs 
with GABAergic interneurons (Ascoli et al. 2008). 

Another possibility is that spatial specificity in the neocortex 
is scale-dependent. It is well known that, at the macroscopic 
and mesoscopic scales, the mammalian nervous system is a 
highly ordered and stereotyped structure where connections 
are established in a highly specific and ordered way. Even at the 
microscopic level, it is also clear that different areas and layers 
of the cortex receive specific inputs. At the ultrastructural level, 
however, our results seem to indicate that synapses are distribu- 
ted in a nearly random pattern. This could mean that, as the 
axon terminals reach their destination, the spatial resolution 
achieved by them is fine enough to find a specific cortical layer, 
but not sufficiently fine to make a synapse on a smaller target 
such as a specific dendritic branch or dendritic spine within 
that layer. For example, axon terminals from thalamic nuclei 
reach specific areas and layers of the cerebral cortex but, once 
there, they would form synapses randomly among their poss- 
ible targets. In other words, at this ultrastructural scale, the 
specificity of connections would not only rely on spatial cues 
but also on other mechanisms such as molecular or activity- 
dependent cues in order to assure synaptic specificity. 

Supplementary Material 

Supplementary material can be found at: http://www.cercor. 
oxfordjournals.org/. 
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